---
title: 'Supporting information: "Consistency in motion event encoding across languages"'
output:
  html_document:
    depth: 2
    number_sections: yes
    theme: default
    toc: yes
    code_folding: hide

---

Set up workspace
================

Libraries:

```{r setup, warning=FALSE, message=FALSE}
library("knitr")
library("tidyverse")
library("scales")     # for y-axis scales as percentages
library("gridExtra")  # for multiplots using arrangeGrob
library("grid")       # grid.draw()
library("brms")       # Bayesian mixed models
library("sjPlot")
library("lme4")
# knitr::opts_chunk$set(echo = TRUE)
```

Load data:

```{r, message=FALSE}
# file paths assume data and script are in the same folder
items <- read_csv("items_info.csv")
descr <- read_csv("descriptions.csv") %>%
  mutate(manner_num = ifelse(manner_expressed == "yes", 1, 0))
```

First and last rows of data file:

```{r}
kable(head(descr, 3))
kable(tail(descr, 3))
```


Source custom functions and set global parameters used in this report:

```{r}
# parameters for ggplots
source("my_ggplot_theme.R")
# functions to plot distributions of DVs by language and speaker/item
source("plot_speakers-events_fnc.R")
# functions to plot entropy in framing choices
source("plot_entropy_fnc.R")
# global parameters
mylangs <- c("Spanish", "Swedish")
myvars  <- c("speaker", "event")
```


Description of stimuli
======================

A full description of the target events used in the study:

```{r}
kable(items)
```


Example calculation of entropy
=============================

The entropy, $H$, of a random variable $X$ is a positive real number which measures that variable's degree of randomness or variability (Cover and Thomas, 2005).
For discrete outcomes as the ones considered here (e.g., the distribution of framing strategies used by a speaker), the entropy of $X$ is defined as 

$$H(X) = - \sum_{x \in X} p(x) \log_2 p(x)$$

where $x$ is any of the possible values that $X$ can take (in our case, 'S', 'V' or 'no path'), and $p(x)$ is its probability.

In the paper, entropy is computed by speakers and by events; the two are computed in a similar way.
We exemplify with by-speaker entropy. 
When computing by-speaker entropy over framing strategies, we are considering each speaker as a random variable.
Let $i$ be a participant in the study.
The variable $X_i$ yields a certain outcome for each description provided by speaker $i$.
The possible outcomes are: V-framing (V), S-framing (S) or no path (nP).
Suppose speaker $i$ described the 32 events using 20 V-descriptions, 10 S-descriptions and 2 descriptions without path; then the entropy $H(X_i)$ of that speaker is estimated as follows:

\begin{equation}
\begin{aligned}
\hat{H}(X_i) &= - \sum_{x \in X} \hat{p}(x) \log_2 \hat{p}(x) \\
&= - \left[ \hat{p}(\text{V}) \log_2 \hat{p}(\text{V}) + \hat{p}(\text{S}) \log_2 \hat{p}(\text{S}) + \hat{p}(\text{nP}) \log_2 \hat{p}(\text{nP}) \right] \\
&= - \left[ \frac{20}{32} \log_2 \frac{20}{32} + \frac{10}{32} \log_2 \frac{10}{32} + \frac{2}{32} \log_2 \frac{2}{32} \right] \\
&= - \left( -0.424 -0.524 -0.250 \right) \\
&= 1.198
\end{aligned}
\end{equation}

where the hat signs in $\hat{H}$ and $\hat{p}$ indicate that we are computing an *estimate* of speaker $i$'s entropy ($\hat{H}$ is sometimes called empirical entropy).
This estimate is based on the maximum-likelihood estimates (here relative frequencies) of that speaker's probability of producing the different outcomes (hence $\hat{p}$).

The entropy of an event $j$ is estimated in an analogous fashion.
Suppose event $j$ was described zero times with a V-description, 31 times with an S-description (S) and once with no path (nP).
The entropy estimate for that event is:

\begin{equation}
\begin{aligned}
\hat{H}(X_j) &= - \sum_{x \in X} \hat{p}(x) \log_2 \hat{p}(x) \\
&= - \left[ \hat{p}(\text{V}) \log_2 \hat{p}(\text{V}) + \hat{p}(\text{S}) \log_2 \hat{p}(\text{S}) + \hat{p}(\text{nP}) \log_2 \hat{p}(\text{nP}) \right] \\
&= - \left[ \frac{0}{32} \log_2 \frac{0}{32} + \frac{31}{32} \log_2 \frac{31}{32} + \frac{1}{32} \log_2 \frac{1}{32} \right] \\
&= - \left( 0 -0.044 -0.156 \right) \\
&= 0.2
\end{aligned}
\end{equation}

Remarks:

- $0 * \log_20$ is defined to be 0 because of the limit calculation
$\lim_{x\to 0} x * log_2(x) = 0$.
- When entropy is computed (as here) by taking the logarithm with base 2, the unit of entropy becomes the bit (other bases are possible, such as taking the natural logarithm).
In this context, bits have the intuitive interpretation of approximately corresponding to the average number of yes/no questions that are required to guess the outcome of the random variable on a specific trial (see Cover and Thomas, 2005, ch. 5).


Code to replicate the reported analyses and figures
============================================

Variability by speakers and items (Fig 4)
-----------------------------

```{r}
# all the work is done by the functions defined in "plot_speakers-events_fnc.R"
mpl_fr <- plot_multiplot(myvars, mylangs, "framing")
grid.draw(mpl_fr)
```


Quantifying variability: entropy analysis
-----------------------------------------

Figure 5:

```{r, fig.height=3}
# see functions defined in "plot_entropy_fnc.R"
set.seed(3122214)
mplot_H_fr <- multiplot_H(descr, "framing")
```

Non-parametric (Mann-Whitney-Wilcoxon) tests comparing entropy by language

```{r}
test_mannwhit(descr, "framing")
```


Variability in manner encoding
------------------------------

Figure 6:

```{r}
mpl_ma <- plot_multiplot(myvars, mylangs, "manner_expressed")
grid.draw(mpl_ma)
```


Expression of manner as a function of framing strategy
----------------------------------------------------

### Figure 7

```{r}
# summarize data
descr_summary <- descr %>%
  group_by(language, framing, manner_expressed) %>%
  summarise(count = n()) %>%
  left_join(descr %>% group_by(language) %>% count()) %>%
  mutate(percent = 100 * count / n)
descr_summary$framing <- factor(descr_summary$framing, levels =  c("V", "S", "No path"))
```

Figure 7:

```{r, fig.height=3}
# plot
ggplot(descr_summary, aes(x = framing, y = percent, fill = manner_expressed)) +
  geom_bar(stat = "identity") +
  ylab("% of descriptions") +
  facet_grid(. ~ language) +
  mytheme
```

As a table (not included in the paper):

```{r}
kable(descr_summary, digits = 1)
```


### Bayesian logistic mixed model of Spanish data

Subset data (Spanish only and exclude 'no path' descriptions):

```{r}
descr_sp <- descr %>% filter(framing != "No path" & language == "Spanish")
```

Define coding scheme:

```{r}
# contrast coding
descr_sp$framing <- factor(descr_sp$framing)
contrasts(descr_sp$framing) <- contr.sum(2)
colnames(contrasts(descr_sp$framing)) <- "S_V"
contrasts(descr_sp$framing)
```

We use the `brm` function from the *brms* package, see e.g. 
[this tutorial](see https://www.rensvandeschoot.com/tutorials/generalised-linear-models-with-brms/).

Fit the model (takes several minutes to run on a laptop):

```{r}
# Comment out either of the following two blocks depending on whether you want
# to run the model from scratch or read it from disk:

# # 1)
# # fit model
# fm_manner_brms <- brm(
#   manner_num ~ 1 + framing +
#     (1 + framing | speaker) + (1 + framing | event),
#   family = bernoulli(link = "logit"),
#   data = descr_sp)
# # write to disk
# write_rds(fm_manner_brms, "fm_manner_brms.rds")

# 2)
# read from disk
fm_manner_brms <- read_rds("fm_manner_brms.rds")
```

Model summary:

```{r}
summary(fm_manner_brms)
```

Brief summary showing population-level (fixed) effects only:

```{r}
tab_model(fm_manner_brms)
```

Plot posterior distributions of population-level (fixed) and random effects:

```{r}
plot(fm_manner_brms)
```

The same but with the same x-axis to enable comparison between effect sizes:

```{r, warning=FALSE}
stanplot(fm_manner_brms, type = "areas", prob = 0.95)
```


Effect of trial on syntactic framing in Spanish descriptions
----------------------------------------------------

A reviewer raised the following issue:

> ... could it be that Spanish speakers produced almost exclusively either
V-descriptions or S-pattern just because they were picking up information as 
they were going through the trials, which would indicate some sort of 
experimental demand characteristic [...]. Perhaps it would be useful to know if
the choice of V-descriptions over S-pattern (or vice versa) was constant over
the course of the experiment for each participant [...].

We address this in a follow-up analysis by fitting a mixed logistic regression
model to the Spanish data.
We model the choice of a verb-framed construction (in log-odds) as a function of
trial.
The model additionally specifies random intercepts by speakers and events and
a random slope for trial by speakers, resulting in the following model 
specification:

`VerbFraming ~ 1 + trial + (1 + trial | speaker) + (1 + trial | event)`


```{r}
descr_sp <- descr_sp %>%
  mutate(
    VF = ifelse(framing == "V", 1, 0),
    trial_c = trial_nb -  mean(trial_nb)  # centre trial
    )
```


```{r}
fm_trial <- glmer(
  VF ~ 1 + trial_c + (1 + trial_c | speaker) + (1 + trial_c | event),
  data = descr_sp,
  family = "binomial"
  )
summary(fm_trial)
```

The model indicates that trial had a numerically positive but non-significant
effect on the log-odds of producing a V-framed description.

We can visually inspect the effect of trial by plotting the best-fitting 
logistic regression curves for the choice of V-pattern as a function of trial.
The following figure fits one regression to the whole data set (thick black
line) and overlays the individual regression lines for each participant 
in light gray:

```{r, message=FALSE}
# plot
ggplot(descr_sp, aes(x = trial_nb, y = VF)) +
  # individual speakers' logistic regression lines
  geom_smooth(aes(group = speaker),
              method = "glm",
              method.args = list(family = "binomial"),
              se = FALSE,
              color = "light gray",
              size = 0.5) +
  # mean logistic regression line
  geom_smooth(method = "glm",
              color = "black",
              method.args = list(family = "binomial")) +
  ylab("Probability of V-framing") +
  xlab("Trial number") +
  mytheme
```

The mean curve has a slight upward slope but the individual lines per participant
show a great deal of variability.
This variability is consistent with the between-speaker variability observed
more generally in the data.

In conclusion, the idea that experimental demands had a consistent effect on
Spanish descriptions is not strongly supported.


Session info
============

```{r}
sessionInfo()
```

